/*    Copyright 2010 Tobias Marschall
 *
 *    This file is part of MoSDi.
 *
 *    MoSDi is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 3 of the License, or
 *    (at your option) any later version.
 *
 *    MoSDi is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoSDi.  If not, see <http://www.gnu.org/licenses/>.
 */

package mosdi.fa;

import java.util.Arrays;

import mosdi.util.BitArray;

public abstract class FiniteMemoryTextModel implements TextModel {
	
	/** Returns the number of states. */
	public abstract int getStateCount();
	
	/** Returns all targets with probability > 0.0. */
	public abstract int[] getTransitionTargets(int sourceState, int character);
	
	/** Probability of going from state sourceState to targetState and thereby
	 *  emitting the given character. */
	public abstract double getProbability(int sourceState, int character, int targetState);
	
	public abstract double getEquilibriumProbability(int state);
	
	public abstract FiniteMemoryTextModel reverseTextModel();
	
	/** Derived classes that do not have an "order", are allowed to return -1. */
	public abstract int order();
	
	/* ****************************************************************************** */
	
	public double[] getEquilibriumDistribution() {
		double[] distribution = new double[getStateCount()];
		for (int state=0; state<getStateCount(); ++state) {
			distribution[state] = getEquilibriumProbability(state);
		}
		return distribution;
	}
	
	/** Probability that, in the given state, the given character is generated. */
	public double getProbability(int sourceState, int character) {
		double p = 0.0;
		for (int targetState : getTransitionTargets(sourceState, character)) {
			p+=getProbability(sourceState, character, targetState);
		}
		return p;
	}

	
	/** Expectation of observing a single occurrence of a pattern assuming to 
	 *  be in equilibrium. */	
	public double expectation(int[] pattern) {
		double[] stateDistribution = stateWiseExpectation(pattern);
		double expectation = 0.0d;
		for (int state = 0; state<getStateCount(); ++state) {
			expectation += stateDistribution[state];
		}
		return expectation;
	}

	/** Assuming to be in equilibrium, returns the expectation of observing a 
	 *  single occurrence of a pattern and ending this occurrence in a specific 
	 *  state. That means, result[c] gives the probability that, starting in 
	 *  equilibrium, the given pattern is found and after the occurrence
	 *  the text model is in state c. */	
	public double[] stateWiseExpectation(int[] pattern) {
		return statewiseProductionProbability(getEquilibriumDistribution(), pattern);
	}

	/** Returns the probability that, starting from initialState, the given
	 *  pattern is produced. */
	public double productionProbability(int initialState, int[] pattern) {
		double result = 0.0;
		for (double d : statewiseProductionProbability(initialState,pattern)) {
			result += d;
		}
		return result;
	}
	
	
	/** Result[state] is the probability that, starting from initialState, the given
	 *  pattern is produced and we end up being in state "state". */
	public double[] statewiseProductionProbability(int initialState, int[] pattern) {
		double[] initialDistribution = new double[getStateCount()];
		initialDistribution[initialState] = 1;
		return statewiseProductionProbability(initialDistribution,pattern);
	}
	
	/** Result[state] is the probability that, starting from initialDistribution, the given
	 *  pattern is produced and we end up being in state "state". */
	public double[] statewiseProductionProbability(double[] initialDistribution, int[] pattern) {
		double[] d0 = initialDistribution;
		for (int c : pattern) {
			double[] d1 = new double[getStateCount()];
			for (int sourceState=0; sourceState<getStateCount(); ++sourceState) {
				if (d0[sourceState]==0.0) continue;
				for (int targetState : getTransitionTargets(sourceState, c)) {
					d1[targetState] += d0[sourceState]*getProbability(sourceState,c,targetState);
				}
			}
			d0 = d1;
		}
		return d0;
	}
	
	/** Expectation of observing a single occurrence of a pattern assuming to 
	 *  be in equilibrium. */	
	public double expectation(BitArray[] pattern) {
		double[] stateDistribution = getEquilibriumDistribution();
		for (BitArray characterSet : pattern) {
			double[] newStateDistribution = new double[getStateCount()];
			if (characterSet.size()!=getAlphabetSize()) throw new IllegalArgumentException("Alphabet size mismatch.");
			for (int c = 0; c < characterSet.size(); ++c) {
				if (!characterSet.get(c)) continue;
				for (int sourceState = 0; sourceState < getStateCount(); ++sourceState) {
					if (stateDistribution[sourceState] == 0.0d) continue;
					for (int targetState : getTransitionTargets(sourceState, c)) {
						newStateDistribution[targetState] += stateDistribution[sourceState] * getProbability(sourceState, c, targetState);	
					}
				}
			}
			stateDistribution = newStateDistribution;
		}
		double expectation = 0.0d;
		for (int state = 0; state<getStateCount(); ++state) {
			expectation += stateDistribution[state];
		}
		return expectation;
	}

	/** Returns, for each of the given generalized characters, the minimal
	 *  possible probability of observing this character. That means, the 
	 *  minimum is taken over all possible states of the text model. 
	 */
	public double[] minProbabilityTable(BitArray[] generalizedAlphabet) {		
		return buildMinOrMaxProbabilityTable(generalizedAlphabet, false);
	}

	/** Returns, for each of the given generalized characters, the maximal
	 *  possible probability of observing this character. That means, the 
	 *  maximum is taken over all possible states of the text model. 
	 */
	public double[] maxProbabilityTable(BitArray[] generalizedAlphabet) {
		return buildMinOrMaxProbabilityTable(generalizedAlphabet, true);
	}
	
	/** Common implementation of buildMinProbabilityTable and buildMaxProbabilityTable. */
	private double[] buildMinOrMaxProbabilityTable(BitArray[] generalizedAlphabet, boolean max) {
		double[] minMaxTable = new double[generalizedAlphabet.length];
		Arrays.fill(minMaxTable, max?0.0d:1.0d);
		for (int i = 0; i < generalizedAlphabet.length; ++i) {        
			for (int state = 0; state < getStateCount(); ++state) { 
				double p = 0.0d;
				for (int c = 0; c < generalizedAlphabet[i].size(); ++c) {
					if (!generalizedAlphabet[i].get(c)) continue;
					p += getProbability(state, c); 
				}
				minMaxTable[i] = max?Math.max(minMaxTable[i], p):Math.min(minMaxTable[i], p);
			}
		}
		return minMaxTable;
	}
	
	/** Returns, for each pair (a,b) of the given generalized characters, the minimal
	 *  possible conditional probability P(S_t=a|S_t=b). That means, the minimum is 
	 *  taken over all possible states of the text model. 
	 *  
	 *   @return result[a][b]<=P(a|b) for all text model states.
	 */
	public double[][] minConditionalProbabilityTable(BitArray[] generalizedAlphabet) {		
		return buildMinOrMaxConditionalProbabilityTable(generalizedAlphabet, false);
	}


	/** Returns, for each pair (a,b) of the given generalized characters, the maximal
	 *  possible conditional probability P(S_t=a|S_t=b). That means, the maximum is 
	 *  taken over all possible states of the text model. 
	 *  
	 *   @return result[a][b]>=P(a|b) for all text model states.
	 */
	public double[][] maxConditionalProbabilityTable(BitArray[] generalizedAlphabet) {		
		return buildMinOrMaxConditionalProbabilityTable(generalizedAlphabet, true);
	}

	/** Common implementation of minConditionalProbabilityTable and maxConditionalProbabilityTable. */
	private double[][] buildMinOrMaxConditionalProbabilityTable(BitArray[] generalizedAlphabet, boolean max) {
		double[][] minMaxTable = new double[generalizedAlphabet.length][generalizedAlphabet.length];
		for (int i=0; i<minMaxTable.length; ++i) {
			Arrays.fill(minMaxTable[i], max?0.0d:1.0d);	
		}
		for (int state = 0; state < getStateCount(); ++state) {
			double[] charProbs = new double[getAlphabetSize()];
			for (int c=0; c<getAlphabetSize(); ++c) charProbs[c] = getProbability(state, c); 
			for (int a = 0; a < generalizedAlphabet.length; ++a) {        
				for (int b = 0; b < generalizedAlphabet.length; ++b) {
					double p1 = 0.0d;
					double p2 = 0.0d;
					for (int c = 0; c < getAlphabetSize(); ++c) {
						if (generalizedAlphabet[b].get(c)) {
							p2+=charProbs[c];
							if (generalizedAlphabet[a].get(c)) p1+=charProbs[c];
						}
					}
					minMaxTable[a][b] = max?Math.max(minMaxTable[a][b], p1/p2):Math.min(minMaxTable[a][b], p1/p2);
				}
			}
		}
		return minMaxTable;
	}

}
